FIFTH ORDER RUNGE-KUTTA-NYSTROM METHODS WITH 
COMPLEX COEFFICIENTS* 
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Abstract. We present fifth order Runge-Kutta-Nystrom methods, where we allow the timestep 
coefficients to assume complex values. Among the methods with complex timesteps, we focus on 

Cn ' the ones with the coefficients that have positive real parts. This property makes them suitable for 

problems where a negative coefficient is not acceptable. In addition, the leading order terms in the 

C ^ , error expansion of these methods are purely imaginary, effectively increasing the order of the methods 

^Nj ' by one for real variables. 
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1. Introduction. Splitting methods [5] provide a number of advantages for 
studying the evolution of Hamiltonian systems. Not only are they simple to im- 
plement, but also they can be fine tuned to exploit the structure of the problem 
at hand [T71 [14] . A large class of physical problems are described by the separable 
r^ ' Hamiltonian 



< 



c^ 



Hiq,p)^T{p) + V{q), (1.1) 



where the kinetic energy T{p) = ^p^ AI~^p is quadratic in momenta, and the potential 
energy V{q) is only a function of the coordinates. For these problems, Hamilton's 
^ , equations 

a^ 

Cn| : dt ~ dp^' dt ~ dqi ^ ■ ' 

m . 

pf-\ ' lead to the second order differential equation 

. ^_, This equation can be efficiently integrated using Runge-Kutta-Nystrom (RKN) meth- 

^ ; ods [ig. 

H ■ RKN methods are particularly effective for orders higher than 4, since some of 

- - - the terms in the error expansion vanish identically, thanks to the special structure 

in equation ([O]) [18l|8]. Zhu & Qin [22] and Okunbor & Skeel ^ found 5th order 
explicit canonical RKN methods with 5 stages, which is the minimum required by 
the order conditions. The stages in these methods have large (around unity) step size 
coefficients, some of which are negative . Large coefficients can lead to large factors for 
integrations errors 0, while negative coefficients can make the method unacceptable 
for the underlying problem |13| . Starting from a splitting scheme. Chambers [llj 
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showed that these difficuhies can be overcome by aUowing the step size coefficients to 
be complex numbers, and found a third order splitting method with two stages. The 
coefficients for this method have small and positive real parts, leading to small errors. 
An interesting property of this method is that the leading order terms in the error 
expansion have purely imaginary coefficients. This makes the method effectively one 
order higher for real variables (see ref. |6], for other examples of such methods). 

In this paper, we present multiple 5th order RKN methods with both real and 
complex step size coefficients. We put special emphasis on methods with coefficients 
that have small, positive real parts, and purely imaginary leading error terms. 

2. 5th Order Canonical RKN Methods. For a system with Hamiltonian 
(jl.ip . we can define "velocity" as g = M^^p. Then an s-stage RKN method is given 
by [Tmi] 



y> 



^ qn + Cihqn + h'^^ o.ijf{yo)^ i = 1, 2, . . . , s , 



(2-1) 

qn+l = qn+ hqn + ^^ X! ^»/(y«)' 9« + l = 9n + ^ X! ^if^Vi) i 

where f{q) is defined in equation (jl.3p . This method is explicit, that is a step depends 
only on previous steps, if aij = for j > i. An explicit s-stage RKN method without 
redundant stages is canonical if [T51 [TO] 

h^B,{l-c,), l<t<s, (2.2) 

a.ij ^ Bj{ci-Cj), 3<i. (2.3) 

For methods of order 5, the following order conditions must hold [19 1 122 ) : 

ii : ^ B, = 1, ^6 : ^ ^ B,B,{c, - c,) = \, 



6' 



h : J2 ^*^' ^l^ *' ^ X E B^B,cUc^ - c,) = 1, (2.4) 

U ■■ Y B,cl = -, *9 ■ E E BiBjCCjia - cj) = — , 

h ■■ Y. ^'^^ = i *io ■ E E E BrBjBi{c, - c,){c, - q) = ^ • 

i i j<i l<i 

By solving these equations, we can obtain an RKN method of order 5; note, 
however, that these equations by themselves do not give us the error behaviour. 

3. Splitting Scheme. Since we will be using a Hamiltonian splitting scheme 
for implementing our integrators and error analysis, we now briefly review the basics 
of this technique using Lie algebra, in a manner similar to the treatment of Yoshida 

First we write Hamilton's equations for a state w{t) — {p{t),q{t)) as 






(3.1) 
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where H, T, V are Lie operators |13] corresponding to the Haniihonian, the kinetic 
energy and the potential energy, respectively. The solution for w{t) can then be 
written as 

w{t) = e"%(0) , (3.2) 

and the exponential operator is approximated by 

e^H«e^^ = ne"'"V'^^. (3.3) 

i 

The expansion for Z in terms of T and V can be obtained by repeatedly applying 
the Baker-Campbell-Hausdorff (BCH) formula (for BCH expressions with minimum 
number of terms, see [H [H])- The coefficients ai, Pi are chosen such that rZ — rH + 
0{t'^'^^), for a method of order n. Note that the expansion of Z can be considerably 
simplified, since the commutator [V, [V, [T, V]]] vanishes for the systems described by 
equation (|1.3I) . 

The evolution due to the exponential operators on the right hand side of equa- 
tion (|3.3I) is given by simple displacements 

e"'"^(g,p) = (g + «,rp,p) and e^-^^(g,p) = (g,p + /3,r/(g)) . (3.4) 

The relation between the coefficients of an RKN method and a splitting scheme is 
given by [18]: 

Qfj = Q - Ci_i, as+i = 1 - Cs,lii = B^, Ps+i = 0, (3.5) 

(where 1 < i < s, cq = 0), leading to the following procedure for the method \TS\ : 

yo = qn, 

yo = Qn, 

for i = 1,2, . . . ,s 

Vi = yi-i + h{ci - Q_i)2)i_i, 

Vi = TJi-i + hBJ{yi), 
qn+i = Vs + h{l - Cs)ys 

g^+i ^ Vs ■ 

This corresponds to a splitting scheme with the structure 
to which we will refer as RKNA scheme. Another method with the structure 



g-rZ ^ g/36-rVga5rTg/35TV ^airT ^I3ir\ 



(3.7) 



to which we will refer as RKNB scheme, would have similar computational cost. How- 
ever, because of the special structure of the RKN methods, the coefficients for the 
latter scheme would be entirely different [5] . To obtain the coefficients from the same 
order conditions, we rewrite equation (13.71) as 



6"^° « e""^ = gaTTTg/JerVgaerTg^srV ^a2TT ^0itV ^uitT /g g\ 

i.e., as a 6 stage method. Since ay = ai = 0, we have ci = and cq = 1. Hence, the 
number of unknowns is again equal to the number of equations. 
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4. Constructing the Methods. We solved the order conditions, eq. (|2.4p . us- 
ing f solve routine of MAPLE, with complex optional keyword. We started from a 
large number of random initial guesses for Bi and Ci in the rectangle —2 < 9{c{z),3m{z) < 
2. This yielded both of the previously known real solutions for RKNA scheme, along 
with their adjoints; as well as three real solutions for RKNB scheme that seems to have 
not been discovered before. In addition, we found a large number of complex solutions 
for both schemes. 

Once the coefficients are calculated, we repeatedly applied BCH formula to cal- 
culate the leading order error terms as nested commutators of T and V. Because of 
the Jacobi identity and the simplification [V, [V, [T, V]]] = 0, not all commutators are 
independent. Also, to calculate the leading order error of the method, we are only 
interested in terms with up to 6 operators; so, we worked in a Philip Hall basis [S] 
with 2 generators and nilpotency 7, leading to 23 terms. Denoting our generators by 
A"i = T and X2 = V, we constructed a "multiplication table" : 

[Xi, ^2] — X3; [Xi^X^] — X4; [A'i,X4] ~ Xq; [XijXr,] = Xj] 

[Xi,Xq] = Ari2; [Ari,Ar7] ^Xg+Xi^-, [Xi,Xg] — Xie, 

[Xi, ATio] = Xii + X17; [Xi, X12] = A'lg; [Xi, ^13] — XiQ + X20; 

[Xi,A'i4] = — A"!! — Xn; [X2,A'3] = A5; [A2,A'4] — Xt, [A'2,A'6] = A'13; (4.1) 

[A2,A'7] = — Xio; [A2,A'g] = —All + A17; [A"2, A12] — A'20; 

[X2, X13] = — 3X17; [X3,X4] = ATg; [Ar3,X5] = Xio; [A'3,X6] ~ Xie; 

[X3, Xr] = Xir, [X4, X5] = Xn; with [X„ X,] = -[X,,X,] . 

Apart from the commutators given in this table, all commutators of the basis operators 
vanish. Using this table, we wrote a Pythor[j program, using the SymPjo package, to 
obtain the expansion for Z up to and including 6 operator terms. This allowed us to 
validate the methods and calculate the coefficients for the leading order error. 

For RKNA scheme, we found two previously published methods with real coeffi- 
cients. For RKNB scheme, we found three methods with real coefficients that we did 
not find elsewhere in the literature. We give the coefficients for these methods in 
Table mH 

We also found a large number of methods with complex coefficients. Among these, 
ten for RKNA scheme and five for RKNB scheme had coefficients with positive real parts 
and purely imaginary leading order errors. We give the coefficients for two methods 
with smallest error coefficients for each scheme in Table 1?^ 

The leading order errors are given by the coefficients in front of the six-term 
commutators in the expansion of Z. These coefficients, for the methods presented, 
are given in Table 14.31 

5. Numerical Experiments. As a simple validation, we first compare the be- 
haviour of two of these methods (RKNACl and RKNBRl) with other methods from the 
literature, for the gravitational two-body problem. We integrated the equations of mo- 
tion of two equal point masses, on orbit around each other with eccentricity e = 0.2, 
for fifty orbital periods. To follow the orbit, we used a Gragg-Bulirsch-Stoer (GBS) 
integrator 115J. At each timestep, we also calculated the expected position and veloc- 
ity for each of the methods we tested. The difference between the outcomes of GBS 
integrator and the method gives an estimate of the error made, as long as they are not 



^ http : //www ■ python . org] 



http : //code . google ■ com/p/sympy | 
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Table 4.1 
Coefficients for real fifth order splitting schemes. The adjoint methods are also valid and can 
be obtained by reversing the order of the coefficients. The coefficients are not entirely independent 
since J^o^ = J]/3j = 1. 

Cti Pi 



ARl 0.96172990014645096 

-0.09525408032034999 

-0.73942683539212613 

0.62730935078241887 

-0.52506178465602220 

0.77070344943962849 

AR2 0.69883375727545265 

-0.49469565362085154 

0.81641946634957295 

-0.65762956677338285 

-0.057841894299102682 

0.69491389106831146 

BRl 0.54200976680171613 

-0.04060817665564392 

-0.87779698530109766 

0.86474236062251646 

0.51165303453250898 

BR2 0.42637413177222316 

-0.82438794434938248 

-0.63140077574154094 

0.38590710518893978 

1.6435074831297605 

BR3 1.0413749845202060 

-0.61784769849171965 

0.62570540985789957 
-0.63446409452971410 

0.58523139864332822 



0.39682804502722538 
-0.824377563589592 
0.2042028689314904 
1.0021847152077973 

0.22116193442307898 

0.40090379269659899 

0.95997088013405985 

0.0884951581272243 

1.2214390923487315 

-1.6708089233066146 

0.24566294009066009 

1.1433587581365421 

-1.3796706973507000 

-0.019611260781217307 

0.87087215441178844 

0.13938810549292669 

0.15102308452230116 

0.72768821316253478 

-0.26217627934521390 

-0.044211509719803855 

0.23596222045571453 

0.19171427092446728 

0.12696076271851077 

-1.4166626058695677 

-0.62172666654176438 

0.69301448863793809 

1.2079876026916669 

1.0104264183632164 



dominated by the truncation errors arising from limited machine accuracy (~ 10^^^). 
We then calculated 0.25 and 0.75 quantile and took their difference, to get the in- 
terquartile range. This is a robust statistic that gives a measure of the dispersion in 
data. 

In figure 15.11 we plot the interquartile range of position error, for various step 
sizes and different methods. For implementing the method of Chambers 11 and our 
methods, we chose to throw away the imaginary part of the positions and the velocities 
after each step. This destroys the symplecticness of the methods (or rather reduces 
the order for which the methods are symplectic) but leads to good error behaviour 

i- 

The comparison indicates that, for this problem, method of Chambers [11] shows 
fourth order and method ACl shows sixth order behaviour, even though they are for- 
mally third and fifth order, respectively. Analogous results were obtained by Cham- 
bers illj for similar problems. It is interesting to see that the optimized 4th order 
method RKNb6 [5] outperforms a method of higher order, RKNbrl, down to machine 
accuracy level. 
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Table 4.2 
Coefficients for a selection of complex fifth order splitting schemes, with purely imaginary lead- 
ing order error terms. All the methods are skew- symmetric, i.e., reversing the order of the co- 
efficients gives their complex conjugates. This property and ^^ cii = ^^ /3 = 1 can be used to 
calculate the missing coefficients. The adjoints of the methods, obtained by reversing the order of 
the coefficients, are also valid. 



ACl 



AC2 



BCl 



BC2 



me 



3m 



0L\ 
(X2 
OL-i 

Pi 
h 

ai 
a2 
as 

Pi 
Pi 

P2 
P3 

ai 
a2 

Pi 

P2 
P3 
Ctl 

a-2. 



0.087808410045663212 
0.17916539354193987 
0.23302619641239692 
0.17526734338348050 
0.18488007701471166 

0.087634204536037057 
0.18007104463252914 
0.23229475083143381 
0.17526840907207411 
0.18487368019298416 

0.093106790861751605 
0.14578332225686154 
0.26110988688138685 
0.15950063058390336 
0.19085044206705213 
0.10625796854753310 
0.35767992721948460 

0.036062104232982296 
0.26934942679787788 
0.14580813747862993 



0.028523844251341822 

-0.067857083007249973 

-0.097952003128893425 

0.057642040076250593 

-0.19410647329733509 

0.028807372065269351 

-0.068253589313355443 

-0.097060961378624794 

0.057614744130538702 

-0.19412192275724959 

-0.026812950639104607 

0.076033669531385746 

0.10851236434561279 

-0.060127448366782494 

0.20369642527600502 

-0.037213537431233983 

-0.022169204268009056 

0.057072185585748646 

-0.093675141997563700 

0.49930185549019606 



Table 4.3 
Coefficients of the terms in the expansion H — Z, for the methods given in TafcZes OTT] and\4-.2\ 



-3.1x10"^ 
5.0x10-3 
1.5x10^4 

-7.0x10-5 
3.5x10-2 

-^2.5x10"*^ 
i3.0xl0"^ 

-iA.lxlO^^ 



4.5x10-* 

-3.1x10-4 

6.2 xlO"-"^ 

9.7x10-4 

1.9x10-4 

il.ixlQ-^ 

-^7.9x10"^ 

-il. 3x10-5 

i4.9xl0-5 



^17 

5.6x10-3 

-5.7x10-3 

-4.1x10-4 

2.7x10-4 

-3.5x10-2 

-i5.0xl0-^ 

ie.OxlQ-^ 

i8.4xl0-*^ 

-a.2xio-4 



-^lE 



1.8x10-5 

1.5x10-5 

5.7x10-5 

5.2x10-4 

-3.0x10-5 

i8.3xl0-^ 

-i9.4xl0-^ 

i2.2xl0-^ 

i7.7xl0-6 



^20 

-6.3x10-5 

2.5x10-4 

1.6x10-4 

5.2x10-4 

-1.0x10-4 

i4.2xl0-s 

-i4.7xl0-6 

-^4.2x10-*^ 

i5.0xl0-5 



JEAx. 



,1/2 



6.44x10-3 
7.58x10-3 
4.71x10-4 
1.25x10-3 
4.90x10-2 
1.02x10-5 
1.15x10-5 
1.96x10-5 
1.49x10-4 



ARl 
AR2 
BRl 
BR2 
BR3 
ACl 
AC2 
BCl 
BC2 



To make a more comprehensive and meaningful test, we also developed integrators 
for the gravitational A^-body problem, based on various methods. We simulated 
the evolution of a Plummer sphere with 400 equal mass particles (see ref. [I] for a 
procedure for constructing a Plummer sphere). We followed the same procedure to 
estimate the errors and calculated the interquartile range for various methods and 
stepsizes. The dependence of errors on step size is given in figure [521 

Since complex arithmetic requires more operations than real arithmetic, we also 
made a comparison of CPU times for various methods. It was not possible to calculate 
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5t/P 



Fig. 5.1. Comparison of the behaviour of various methods with respect to changing step size, 
for gravitational two-body problem, x-axis the logarithm of the time step measured in orbital periods, 
y-axis is the logarithm of the interquartile range for position error. Leapfrog (down-triangles) is 
the standard second order Stormer-Verlet scheme, Triple jump (pluses) is the fourth order scheme 
of Yoshida J21f . Chambers (up-triangles) is the complex third order scheme of Chambers Jllf .RKHhS 
(stars) and RKNal4 (circles) are optimized fourth and sixth order methods by Blanes & Moan f^ [Sp, 
RKNbrl (pentagons) and RKNacl (squares) are two of the methods presented in this work. 



Table 5.1 
CPU times spent at each step for various methods. 
CPU time. 



Values are normalized to Leapfrog method's 



Leapfrog 


1 


Chcimbers 


12 


Triple Jump 


3 


RKNb6 


6 


RKNal4 


14 


RKNbrl 


5 


RKNacl 


28 



CPU times for two-body or 400-body problems accurately, since each step took too 
little time. Consequently, we setup a system of 10000 particles and integrated over 
a few steps. While calculating the CPU time, we subtracted off the time spent for 
the last substep for an integration step, since all the schemes we consider have the so 
called first-same-as-last property. For example, in RKNA schemes (equation l3.6p e^i"^^ 
substep can be combined with e"**'^'"' substep of the next step. We present the CPU 
times per step for various methods in table 15.11 

The data here show that the integration time is proportional to the number of 
stages and using complex arithmetic increases the computation cost by about a factor 
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Leapfrog 

Chambers 

RKNb6 

RKNal4 

RKNbrl 

RKNacl 



10 



-2 



10- 



St (code units) 



Fig. 5.2. Comparison of the behaviour of various methods with respect to changing step size, 
for gravitational 400-body problem, x-axis the logarithm of the time step in code units, y-axis is the 
logarithm of the interquartile range for position error. The symbols correspond to the same methods 
as fiaure l57l\ except Triple Jump is omitted, since it nearly coincides with Chambers for this problem. 



of 6. This last factor surely depends on the problem, the implementation and the 
compiler. We present the error vs. CPU time, based on these findings in figure [5731 



6. Discussion. In this work we constructed fifth order RKN methods with com- 
plex coefficients. Some of the methods we found have much smaller timesteps than 
the previously known methods, leading to smaller leading order errors. In addition, 
many of them have positive real coefficients, making them suitable for problems where 
negative real timesteps are not acceptable [12l[3j[T3]. Similar methods, satisfying this 
requirement, were also developed by Bandrauk et al. [2]. Other high order methods 
with complex coefficients with small positive real parts were developed and analyzed 
by Hansen & Ostermann [16] and Castella et al. [10]; however these authors did 
not specifically study the RKN case, which allows considerable simplification for high 
order schemes. 

For problems where complex arithmetic and/or positive real parts of the coeffi- 
cients is a necessity, we expect these methods to be already competitive with lower 
order methods. However, the results of Blanes & Moan [8 and Sofroniou & Spaletta 
[20| suggest that there is room for improvement by increasing the number of stages. 
We consider finding fully optimized methods with more stages beyond the scope of 
this paper; such an effort would require developing new software and would need 
spending considerable computational power. However, the methods found here can 
be improved by readily available tools. The idea is to turn a 5-stage method into a 
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pa 
o 

I 

PI 
03 



IQ- 



10- 



10 



10- 



-6 



CD 



IQ- 



IQ- 



10 



-14 




Leapfrog 

Triple Jump 

Chambers 

RKNb6 

RKNal4 

RKNbrl 

RKNacl 



10-1 10" 10^ 

CPU time (normalized to LF, St = 0.01) 



10^ 



Fig. 5.3. CPU time spent vs. error for various methods applied to gravitational 400-body 
problem, x-axis the logarithm of the CPU time, normalized to time spent for Leapfrog method at 
time .step St = 0.01, y-axis is the logarithm of the interquartile range for position error. The symbols 
correspond to the same methods as figure I5.il 



6-stage one by adding an additional stage. 



,airT P[t\ 



(6.1) 



and set f3'^ = /3[ = 0. We can then start in the vicinity of this solution and use 
maple's minimization routines. Starting from solution ACl in table 14.21 we obtained 
the following skew-symmetric method: 



ai = .101907705405177865 + i 
a2 = .218628781976265590 -l-i 
as = .179463512618556560- i 
P[ = .0489489561074426954 + 
^2 = .166479171860817010+ i 
^^ = .192297943665939275- 



.130701756906677735, 
.0126440811480678494, 
.148112326926992222, 
i.0669384556781967844, 
.0764027877516731402 , 
.0835834606213808479 , 



(6.2) 



f3'. 



.184547856731601789. 



,1/2 



This method has an error {J2- \Xi\^)'''^ = 2.6 x 10 ^, about two orders of magnitude 



smaller than the other methods (see table 14.31) . 

Since we minimized the (imaginary) sixth order error terms, this optimization 
does not benefit the solutions of gravitational A^-body problem. A possible venue 
for exploration is to increase the number of stages and minimize the (real part of) 
seventh order errors by using the extra variables. For this investigation, using a 
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Lyndon basis would be more preferable, since the BCH expansion has much fewer 
terms in this basis [5] and the consequences of the simplification [V, [V, [T, V]]] = 
is more straightforward. Here, we used a Philip Hall basis, since this allowed us to 
check our algebra by the "Lie Tools Package"|f| of Miguel Torres- Torriti. 

The source code of the MAPLE, Python and C programs used in this work are freely 



available online: https : //github/ataLkaxL/Complex_Coef f _RKN| 
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